{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/glade/work/tjuliano/conda-envs/myclone/lib/python3.9/site-packages/pandas/core/computation/expressions.py:21: UserWarning: Pandas requires version '2.8.4' or newer of 'numexpr' (version '2.8.3' currently installed).\n",
      "  from pandas.core.computation.check import NUMEXPR_INSTALLED\n"
     ]
    }
   ],
   "source": [
    "from netCDF4 import Dataset\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from matplotlib.cm import get_cmap\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import matplotlib.gridspec as gridspec\n",
    "from wrf import getvar, interplevel, to_np, get_basemap, latlon_coords, destagger\n",
    "import glob\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "import cmocean\n",
    "import pyresample\n",
    "\n",
    "### Import the scripts needed for importing satellite data ###\n",
    "from functions_import import import_VNP02, import_cldprop_l2, import_cldprop_uc_l2, import_MOD021KM, import_MOD06_L2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = \"mynn\"\n",
    "dom1 = \"d01\"\n",
    "dom2 = \"d02\"\n",
    "active_dom = \"d02\"\n",
    "sens = \"no_fractional_seaice\"\n",
    "day = \"12\"\n",
    "date = \"0313\"\n",
    "#root = \"/glade/scratch/tjuliano/doe_comble/new_lam_runs/WRF/WRF-ASR-CAO/test/\"\n",
    "root = \"/glade/campaign/uwyo/wyom0122/lam/new_runs/\"\n",
    "cases = [\"031220_no_edmf\",\"031220_mynn_l3\",\"031220\",\"031220_ysu\"]\n",
    "rst = \"/RESTART_D\"\n",
    "if active_dom == \"d01\":\n",
    "    plt_titles = [\"Level 2.5\", \"Level 3\", \"Level 2.5 EDMF\", \"YSU\"]\n",
    "    plt_titles2 = [\"Level 2.5\", \"Level 3\", \"Level 2.5 EDMF\", \"YSU\"]\n",
    "    dx = 3000.\n",
    "else:\n",
    "    plt_titles = [\"Level 2.5\", \"Level 3\", \"Level 2.5 EDMF\", \"YSU\"]\n",
    "    plt_titles2 = [\"Level 2.5\", \"Level 3\", \"Level 2.5 EDMF\", \"YSU\"]\n",
    "    dx = 1000.\n",
    "savedir = \"figs_new_nz136/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lwprng = np.arange(-30,80,5)\n",
    "times = [9,11,13,19,25,31,37,43,49]\n",
    "time_str = ['28_1030z','28_1130z','28_1230z','28_1530z','28_1830z','28_2130z','29_0030z','29_0330z','29_0630z']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "lev = 10\n",
    "skipx = 6\n",
    "skipy = 6\n",
    "skipx2 = skipx*3\n",
    "skipy2 = skipy*3\n",
    "scale = 50\n",
    "barbw = 0.04\n",
    "g = 9.81"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = Dataset(root + cases[0] + \"/wrfinput_\" + dom1)\n",
    "cosalpha = s1.variables['COSALPHA'][-1,:,:]\n",
    "sinalpha = s1.variables['SINALPHA'][-1,:,:]\n",
    "s1.close()\n",
    "\n",
    "s1 = Dataset(root + cases[3] + \"/wrfinput_\" + dom2)\n",
    "cosalpha2 = s1.variables['COSALPHA'][-1,:,:]\n",
    "sinalpha2 = s1.variables['SINALPHA'][-1,:,:]\n",
    "s1.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "region = 'cells'\n",
    "if region == 'cells':\n",
    "###### OPEN CELLS\n",
    "#\tlon_arr = [12.0, 12.0, 16.0, 16.0]\n",
    "#\tlat_arr = [69.5, 71.5, 71.5, 69.5]\n",
    "\tlon_arr = [12.5, 12.5, 16.5, 16.5]\n",
    "\tlat_arr = [69.0, 71.0, 71.0, 69.0]\n",
    "elif region == 'rolls':\n",
    "###### ROLLS/CLOSED CELLS\n",
    "##        lon_arr = [7.5,  7.5, 12.5, 12.5]\n",
    "##        lat_arr = [71.5, 73.0, 73.0, 71.5]\n",
    "\tlon_arr = [8.0, 8.0, 12.0, 12.0]\n",
    "\tlat_arr = [74.5, 76.5, 76.5, 74.5]\n",
    "else:\n",
    "        print ('You are an idiot!')\n",
    "        \n",
    "#lat2 = getvar(s1,\"lat\",meta=False)\n",
    "#lon2 = getvar(s1,\"lon\",meta=False)\n",
    "#hgt = getvar(s1,\"HGT\",meta=False)\n",
    "#lat_d01 = lat2\n",
    "#lon_d01 = lon2\n",
    "#hgt_d01 = hgt\n",
    "ll_lon = lon_arr[0]\n",
    "ll_lat = lat_arr[0]\n",
    "ul_lon = lon_arr[0]\n",
    "ul_lat = lat_arr[1]\n",
    "ur_lon = lon_arr[2]\n",
    "ur_lat = lat_arr[1]\n",
    "lr_lon = lon_arr[2]\n",
    "lr_lat = lat_arr[0]\n",
    "# center point of overall domain\n",
    "ref_lat = (ll_lat + ur_lat)/2.\n",
    "#ref_lon = -(abs(ll_lon) + abs(ur_lon))/2.\n",
    "ref_lon = (ll_lon + ur_lon)/2.\n",
    "######################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Doing  ['VJ102IMG', 'VJ103IMG', 'CLDPROP_L2_VIIRS_NOAA20', 'VJ103MOD', 'viirs_temp']\n",
      "Doing  1212\n",
      "/glade/scratch/tjuliano/doe_comble/satellite/data/0313/VJ102IMG/VJ102IMG\n"
     ]
    }
   ],
   "source": [
    "### Head directory\n",
    "head = \"/glade/scratch/tjuliano/doe_comble/satellite/data/\" + date + \"/\"\n",
    "\n",
    "### Satellite file names, --> Imagery, Imagery geo location, cloud retrival, cloud retrival geo location, folder name ###\n",
    "sats = [[\"VJ102IMG\",\"VJ103IMG\",\"CLDPROP_L2_VIIRS_NOAA20\", \"VJ103MOD\", \"viirs_temp\"]] # NOAA-20\n",
    "\n",
    "### Times of the satellite overpasses in UTC###\n",
    "times = [[\"1212\"]] # NOAA-20\n",
    "\n",
    "### Corresponding WRF times in UTC###\n",
    "wrf_times = [[\"12:00\"]] # NOAA-20\n",
    "\n",
    "wrf_idx = [[12]]\n",
    "\n",
    "### Latitudes of domain center ###\n",
    "lats = [[71.70]] # NOAA-20\n",
    "\n",
    "### Longitudes of domain center ###\n",
    "lons = [[12.80]] # NOAA-20\n",
    "\n",
    "###  Loop over all satellites ###\n",
    "sat_i = -1 # index to get satellite\n",
    "\n",
    "for sat in sats:\n",
    "    print ('Doing ', sat)\n",
    "    sat_i += 1 # increase to get next satellite\n",
    "    coord_i = -1 # Index to get coordinate\n",
    "\n",
    "    count2 = 0\n",
    "    ### Loop over all times ###\n",
    "    for time in times[sat_i]:\n",
    "        print ('Doing ', time)\n",
    "            \n",
    "        coord_i += 1 # increase coordinate index\n",
    "        lon_0, lat_0 = lons[sat_i][coord_i], lats[sat_i][coord_i] # lon and lat at this time\n",
    "\n",
    "        ### all files ###\n",
    "        file = sorted(glob.glob(head+sat[0]+\"/\"+sat[0]+\".A2020073.\"+time+\"*\"))\n",
    "        print (head+sat[0]+\"/\"+sat[0])\n",
    "        geo_file = sorted(glob.glob(head+sat[1]+\"/\"+sat[1]+\".A2020073.\"+time+\"*\"))\n",
    "        cloud_file = sorted(glob.glob(head+sat[2]+\"/\"+sat[2]+\".A2020073.\"+time+\"*\"))\n",
    "        cloud_geo_file = sorted(glob.glob(head+sat[3]+\"/\"+sat[3]+\".A2020073.\"+time+\"*\"))\n",
    "\n",
    "        ### Import sat files using external script ###\n",
    "        if sat[0][0] == \"V\": # VIIRS\n",
    "            data, lon, lat = import_VNP02.import_VNP02(file[0],geo_file[0]) # Imagery\n",
    "            cth, ctt_sat, cwp, cot, cer, cph, cloud_lon, cloud_lat = import_cldprop_l2.import_cldprop_l2(cloud_file[0], cloud_geo_file[0]) # cloud retrivals\n",
    "            cth_uc, ctt_uc, cwp_uc, cot_uc, cer_uc, cloud_uc_lon, cloud_uc_lat = import_cldprop_uc_l2.import_cldprop_uc_l2(cloud_file[0], cloud_geo_file[0]) # cloud retrivals\n",
    "        if sat[0][0] == \"M\": # MODIS\n",
    "            data, lon, lat = import_MOD021KM.import_MOD021KM(file[0],geo_file[0]) # Imagery\n",
    "            cth, ctt_sat, cwp, cot, cer, cph, cloud_lon, cloud_lat = import_MOD06_L2.import_MOD06_L2(cloud_file[0], cloud_geo_file[0]) # cloud retrivals\n",
    "            \n",
    "        ### Calculate coordinates of corner of LES domain ###\n",
    "        #g = pyproj.Geod(ellps='WGS84')\n",
    "        #diagonal = (2*(frame_size**2))**(0.5)/2 # distance from midpoint of domain to corners\n",
    "        #urc = np.asarray(g.fwd(lon_0, lat_0, 45, diagonal)[0:2]) # upper right corner [lon, lat]\n",
    "        #lrc = np.asarray(g.fwd(lon_0, lat_0, 135, diagonal)[0:2]) # lower right corner [lon, lat]\n",
    "        #llc = np.asarray(g.fwd(lon_0, lat_0, 225, diagonal)[0:2]) # lower left corner [lon, lat]\n",
    "        #ulc = np.asarray(g.fwd(lon_0, lat_0, 315, diagonal)[0:2]) # upper left corner [lon, lat]\n",
    "\n",
    "        ### Mask cloud retrival pixels outside LES domain ###\n",
    "        ### This cloud mask is not perfect especially if the domain is larger ###\n",
    "        ### Because the longitude changes along the sides ###\n",
    "        #mask_cloud = np.logical_or.reduce((cloud_lon>(urc[0]+lrc[0])/2, cloud_lat>urc[1], cloud_lon<(ulc[0]+llc[0])/2, cloud_lat<llc[1]))\n",
    "        mask_cloud = np.logical_or.reduce((cloud_lon>(ur_lon+lr_lon)/2, cloud_lat>ur_lat, cloud_lon<(ul_lon+ll_lon)/2, cloud_lat<ll_lat))\n",
    "        cth[mask_cloud],ctt_sat[mask_cloud],cwp[mask_cloud],cer[mask_cloud],cot[mask_cloud],cph[mask_cloud] = np.nan, np.nan, np.nan, np.nan, np.nan, np.nan\n",
    "        cth_uc[mask_cloud],ctt_uc[mask_cloud],cwp_uc[mask_cloud],cer_uc[mask_cloud],cot_uc[mask_cloud] = np.nan, np.nan, np.nan, np.nan, np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading file...\n",
      "/glade/campaign/uwyo/wyom0122/lam/new_runs/031220_no_edmf/RESTART_D/wrfout_wind_d02_2020-03-13_12_00_00\n",
      "[104 326 342 121] [722 705 849 881]\n",
      "Interpolating...\n",
      "Plotting LWP...\n",
      "Reading file...\n",
      "/glade/campaign/uwyo/wyom0122/lam/new_runs/031220_mynn_l3/RESTART_D/wrfout_wind_d02_2020-03-13_12_00_00\n",
      "Interpolating...\n",
      "Plotting LWP...\n",
      "Reading file...\n",
      "/glade/campaign/uwyo/wyom0122/lam/new_runs/031220/RESTART_D/wrfout_wind_d02_2020-03-13_12_00_00\n",
      "Interpolating...\n",
      "Plotting LWP...\n",
      "Reading file...\n",
      "/glade/campaign/uwyo/wyom0122/lam/new_runs/031220_ysu/RESTART_D/wrfout_wind_d02_2020-03-13_12_00_00\n",
      "Interpolating...\n",
      "Plotting LWP...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/glade/work/tjuliano/casper_npl_clone/lib/python3.6/site-packages/ipykernel_launcher.py:175: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.\n"
     ]
    }
   ],
   "source": [
    "files_wind0 = sorted(glob.glob(root + cases[0] + rst + \"/wrfout_wind_\" + active_dom + \"*\"))\n",
    "files_wind1 = sorted(glob.glob(root + cases[1] + rst + \"/wrfout_wind_\" + active_dom + \"*\"))\n",
    "files_wind2 = sorted(glob.glob(root + cases[2] + rst + \"/wrfout_wind_\" + active_dom + \"*\"))\n",
    "files_wind3 = sorted(glob.glob(root + cases[3] + rst + \"/wrfout_wind_\" + active_dom + \"*\"))\n",
    "\n",
    "files_cloud0 = sorted(glob.glob(root + cases[0] + rst + \"/wrfout_cloud_\" + active_dom + \"*\"))\n",
    "files_cloud1 = sorted(glob.glob(root + cases[1] + rst + \"/wrfout_cloud_\" + active_dom + \"*\"))\n",
    "files_cloud2 = sorted(glob.glob(root + cases[2] + rst + \"/wrfout_cloud_\" + active_dom + \"*\"))\n",
    "files_cloud3 = sorted(glob.glob(root + cases[3] + rst + \"/wrfout_cloud_\" + active_dom + \"*\"))\n",
    "\n",
    "optdep_thresh = 0.1\n",
    "method = 'huang'\n",
    "if method == 'huang':\n",
    "    # huang et al. (2014) following Dudhia (1989)\n",
    "    abscoefw = 0.145\n",
    "    abscoefr = 0.00033\n",
    "    abscoefi = 0.0735\n",
    "    abscoefs = 0.00234\n",
    "elif method == 'wrf-python':\n",
    "    # wrf-python\n",
    "    abscoefi = 0.272\n",
    "    abscoefw = 0.145\n",
    "else:\n",
    "    print ('You are stupid!')\n",
    "\n",
    "times = [12]\n",
    "first_pass1 = 0\n",
    "for i in times: #np.arange(12,len(files_cloud0),6):\n",
    "\tfor j in np.arange(len(cases)):\n",
    "\t\tprint ('Reading file...')\n",
    "\t\tif j == 0:\n",
    "\t\t\ts1 = Dataset(files_wind0[i])\n",
    "\t\t\ts2 = Dataset(files_cloud0[i])\n",
    "\t\t\tprint (files_wind0[i])\n",
    "\t\telif j == 1:\n",
    "\t\t\ts1 = Dataset(files_wind1[i])\n",
    "\t\t\ts2 = Dataset(files_cloud1[i])\n",
    "\t\t\tprint (files_wind1[i])\n",
    "\t\telif j == 2:\n",
    "\t\t\ts1 = Dataset(files_wind2[i])\n",
    "\t\t\ts2 = Dataset(files_cloud2[i])\n",
    "\t\t\tprint (files_wind2[i])\n",
    "\t\telif j == 3:\n",
    "\t\t\ts1 = Dataset(files_wind3[i])\n",
    "\t\t\ts2 = Dataset(files_cloud3[i])\n",
    "\t\t\tprint (files_wind3[i])\n",
    "\t\tif j == 0:\n",
    "\t\t\ts4 = Dataset(root + cases[0] + \"/wrfinput_\" + active_dom)\n",
    "\t\t\tlat = getvar(s4,\"lat\",meta=False)\n",
    "\t\t\tlon = getvar(s4,\"lon\",meta=False)\n",
    "\n",
    "\t\tif j == 0:\n",
    "\t\t\tfig = plt.figure(figsize=(18,15))\n",
    "\t\t\tax1 = fig.add_axes([0.7,0.05,0.05,0.4])\n",
    "        \n",
    "\t\tif first_pass1 == 0:\n",
    "\t\t\tlon2d = s1.variables['XLONG'][-1,:,:]\n",
    "\t\t\tlat2d = s1.variables['XLAT'][-1,:,:]\n",
    "\t\t\tgrid = pyresample.geometry.GridDefinition(lats=lat2d, lons=lon2d)\n",
    "\n",
    "            # Define points\n",
    "\t\t\tswath = pyresample.geometry.SwathDefinition(lons=lon_arr, lats=lat_arr)\n",
    "            # Determine nearest (w.r.t. great circle distance) neighbour in the grid.\n",
    "\t\t\t_, _, index_array, distance_array = pyresample.kd_tree.get_neighbour_info(\n",
    "\t\t\tsource_geo_def=grid, target_geo_def=swath, radius_of_influence=5000,neighbours=1)\n",
    "            # get_neighbour_info() returns indices in the flattened lat/lon grid. Compute the 2D grid indices:\n",
    "\t\t\tcalipso_j1, calipso_i1 = np.unravel_index(index_array, grid.shape)\n",
    "\t\t\tprint (calipso_j1, calipso_i1)\n",
    "\n",
    "\t\t\tfirst_pass1 = 1\n",
    "\n",
    "\t\tcalipso_j = calipso_j1\n",
    "\t\tcalipso_i = calipso_i1\n",
    "\n",
    "\t\tw = getvar(s1,\"wa\",units=\"m s-1\",meta=False)\n",
    "\t\ttemp = getvar(s1,\"temp\",units='K',meta=True)\n",
    "\t\ttemp2 = getvar(s1,\"temp\",units='K',meta=False)\n",
    "\t\tpres = getvar(s1,\"p\",units=\"hPa\",meta=True)\n",
    "\t\tpres2 = getvar(s1,\"p\",units=\"Pa\",meta=False)\n",
    "\t\tdp = -1.*(pres2[1:,:,:]-pres2[0:-1,:,:])\n",
    "\t\tdp = np.insert(dp,0,dp[0,:,:],axis=0)\n",
    "\t\tqvapor = s2.variables['QVAPOR'][-1,:,:,:]*1000.\n",
    "\t\tqcloud = s2.variables['QCLOUD'][-1,:,:,:]*1000.\n",
    "\t\tqice = s2.variables['QICE'][-1,:,:,:]*1000.\n",
    "\t\tqsnow = s2.variables['QSNOW'][-1,:,:,:]*1000.\n",
    "\t\tqrain = s2.variables['QRAIN'][-1,:,:,:]*1000.\n",
    "\t\tqcloud2 = s2.variables['QCLOUD'][-1,:,calipso_j[0]:calipso_j[2],calipso_i[1]:calipso_i[3]]*1000.\n",
    "\t\tqice2 = s2.variables['QICE'][-1,:,calipso_j[0]:calipso_j[2],calipso_i[1]:calipso_i[3]]*1000.\n",
    "\t\tqsnow2 = s2.variables['QSNOW'][-1,:,calipso_j[0]:calipso_j[2],calipso_i[1]:calipso_i[3]]*1000.\n",
    "\t\tqrain2 = s2.variables['QRAIN'][-1,:,calipso_j[0]:calipso_j[2],calipso_i[1]:calipso_i[3]]*1000.\n",
    "\t\tph = getvar(s1,\"PH\",meta=True)\n",
    "\t\tphb = getvar(s1,\"PHB\",meta=True)\n",
    "\t\tgph = (ph + phb) / g\n",
    "\t\tgph = destagger(gph,0)\n",
    "\t\tter = getvar(s1,\"HGT\",meta=True)\n",
    "\t\tz = getvar(s1,\"z\", units=\"m\",meta=False)\n",
    "\t\tpblh = getvar(s1,\"PBLH\",meta=False)\n",
    "        \n",
    "\t\tprint ('Interpolating...')\n",
    "\t\thgt = pblh/2.\n",
    "\t\twi = interplevel(w, z, hgt)\n",
    "\n",
    "\t\tliq_cond = qcloud # + qrain\n",
    "\t\tfrz_cond = qice # + qsnow + qgraup\n",
    "\n",
    "\t\tif method == 'huang':\n",
    "\t\t\toptdep = (abscoefw*liq_cond + abscoefr*qrain + abscoefi*qice + abscoefs*qsnow)*dp/g\n",
    "\t\t\toptdep2 = optdep[:,calipso_j[0]:calipso_j[2],calipso_i[1]:calipso_i[3]]\n",
    "                \n",
    "\t\t\ttemp3 = temp2[:,calipso_j[0]:calipso_j[2],calipso_i[1]:calipso_i[3]]\n",
    "                \n",
    "\t\t\tctt2 = np.zeros([np.shape(optdep2)[1],np.shape(optdep2)[2]])\n",
    "\t\t\tfor qq in np.arange(np.shape(optdep2)[1]):\n",
    "\t\t\t\tfor ww in np.arange(np.shape(optdep2)[2]):\n",
    "\t\t\t\t\thold = 0.0\n",
    "\t\t\t\t\tfor ee in np.arange(np.shape(optdep2)[0]-1,0,-1):\n",
    "\t\t\t\t\t\thold = hold + (abscoefw*qcloud2[ee,qq,ww] + abscoefr*qrain2[ee,qq,ww] + abscoefi*qice2[ee,qq,ww] + abscoefs*qsnow2[ee,qq,ww])*dp[ee,qq,ww]/g\n",
    "\t\t\t\t\t\tif hold > optdep_thresh:\n",
    "\t\t\t\t\t\t\tctt2[qq,ww] = temp3[ee,qq,ww] - 273.15\n",
    "\t\t\t\t\t\t\tbreak\n",
    "\t\t\t\t\t\tif ee == 1:\n",
    "\t\t\t\t\t\t\tctt2[qq,ww] = np.nan\n",
    "\t\t\t\t\t\t\tbreak\n",
    "                \n",
    "\t\telif method == 'wrf-python':\n",
    "\t\t\tctt2 = ctt(pres,temp,qvapor,liq_cond,gph,ter,qice=frz_cond,fill_nocloud=True,missing=999.9,opt_thresh=optdep_thresh,meta=False)\n",
    "\t\t\tctt2 = ctt2[calipso_j[0]:calipso_j[2],calipso_i[1]:calipso_i[3]]\n",
    "\n",
    "\t################################################################################################################\n",
    "\t############################################### LWP ############################################################\n",
    "\t################################################################################################################\n",
    "\t\tprint ('Plotting LWP...')\n",
    "        \n",
    "\t\tplt.subplot(2,3,j+1)\n",
    "\t\tm = Basemap(llcrnrlon=ll_lon,llcrnrlat=ll_lat,urcrnrlon=ur_lon,urcrnrlat=ur_lat,lon_0=ref_lon,lat_0=ref_lat,projection='cass',resolution='h')\n",
    "\t\tx, y = m(to_np(lon), to_np(lat))\n",
    "\t\t# Make the contour plot\n",
    "\t\tnorm = mpl.colors.Normalize(vmin=0, vmax=3000)\n",
    "\t\trng = np.arange(-50,-14,2)\n",
    "\t\tctt_contourf = m.contourf(x[calipso_j[0]:calipso_j[2],calipso_i[1]:calipso_i[3]],y[calipso_j[0]:calipso_j[2],calipso_i[1]:calipso_i[3]],to_np(ctt2[:,:]),cmap=cmocean.cm.thermal,levels=rng,extend='both',zorder=2)\n",
    "\t\tw_contours = m.contour(x,y,to_np(wi[:,:]),colors='magenta',levels=[0.25],zorder=3)\n",
    "\n",
    "\t\t# Add geographic outlines\n",
    "\t\tm.drawcoastlines(linewidth=2.5,zorder=3)\n",
    "\t\tif j == 0 or j == 3:\n",
    "\t\t\tm.drawparallels(np.arange(60.,94.,0.5),labels=[True,False,False,False],fontsize=22,linewidth=3.0,zorder=3)\n",
    "\t\t\tm.drawmeridians(np.arange(-40.,70.,1.5),labels=[False,False,False,False],fontsize=22,linewidth=3.0,zorder=3)\n",
    "\t\t\tif j == 3:\n",
    "\t\t\t\tm.drawmeridians(np.arange(-40.,70.,1.5),labels=[False,False,False,True],fontsize=22,linewidth=3.0,zorder=3)\n",
    "\t\telse:\n",
    "\t\t\tm.drawparallels(np.arange(60.,94.,0.5),labels=[False,False,False,False],fontsize=22,linewidth=3.0,zorder=3)\n",
    "\t\t\tm.drawmeridians(np.arange(-40.,70.,1.5),labels=[False,False,False,False],fontsize=22,linewidth=3.0,zorder=3)\n",
    "\t\tm.drawcountries(linewidth=2.0,zorder=3)\n",
    "\t\tm.fillcontinents(color='coral',lake_color='coral',zorder=3)\n",
    "        \n",
    "\t\tx_anx,y_anx = m(15.874211,69.232286)\n",
    "\t\tm.scatter(x_anx,y_anx,marker='D',s=240,facecolor='yellow',edgecolor='k',zorder=8)\n",
    "        \n",
    "\t\tplt.title(plt_titles[j],fontdict = {'fontsize' : 22})\n",
    "\t\tif j == 1:\n",
    "\t\t\tplt.suptitle(files_wind0[i][-19:],y=1.02,fontsize=28)\n",
    "                \n",
    "\t\tif j == 3:\n",
    "\t\t\tplt.tight_layout()\n",
    "\t\t\tplt.subplots_adjust(right=0.90)\n",
    "\t\t\tcbar = fig.colorbar(ctt_contourf, cax=ax1)\n",
    "\t\t\tcbar.ax.tick_params(labelsize=18)\n",
    "\t\t\tcbar.ax.set_ylabel('CTT (K)', labelpad=40, rotation=270, size=24)\n",
    "   \n",
    "\tplt.subplot(2,3,5)\n",
    "\tm = Basemap(llcrnrlon=ll_lon,llcrnrlat=ll_lat,urcrnrlon=ur_lon,urcrnrlat=ur_lat,lon_0=ref_lon,lat_0=ref_lat,projection='cass',resolution='h')\n",
    "\tctt_contourf = m.contourf(cloud_lon, cloud_lat, ctt_sat,latlon=True,cmap=cmocean.cm.thermal,levels=rng,extend='both',zorder=2)\n",
    "\n",
    "\t# Add geographic outlines\n",
    "\tm.drawcoastlines(linewidth=2.5,zorder=3)\n",
    "\tm.drawparallels(np.arange(60.,94.,0.5),labels=[False,False,False,False],fontsize=22,linewidth=3.0,zorder=3)\n",
    "\tm.drawmeridians(np.arange(-40.,70.,1.5),labels=[False,False,False,True],fontsize=22,linewidth=3.0,zorder=3)\n",
    "\tm.drawcountries(linewidth=2.0,zorder=3)\n",
    "\tm.fillcontinents(color='coral',lake_color='coral',zorder=3)\n",
    "        \n",
    "\tx_anx,y_anx = m(15.874211,69.232286)\n",
    "\tm.scatter(x_anx,y_anx,marker='D',s=240,facecolor='yellow',edgecolor='k',zorder=8)\n",
    "        \n",
    "\tplt.title('NOAA-20 VIIRS',fontdict = {'fontsize' : 22})\n",
    "\n",
    "\tplt.savefig(savedir + \"compare_mar\" + day + \"_mynn_edmf_zoom_more_\" + region + \"_\" + active_dom + \"_PANEL_CTT_with_sat_0p5zi_031320_1800.png\",dpi=200,bbox_inches='tight')\n",
    "\tplt.close()\n",
    "    \n",
    "\ts1.close()\n",
    "\ts2.close()\n",
    "\ts4.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
